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Geometric properties of acoustic waves generated by a point source 
in solar-like interior: effects of acoustic cut-off frequency. 
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ABSTRACT 

Acoustic waves generated by a point source in stratified plasma are considered in this paper 
Analytical parametric solution for monochromatic source is derived for plane-parallel poly- 
trope model of the solar interior The solution is used to gain insight into the properties of 
the generated wavefront as a function of excitation frequency and depth. A slowly varying 
pressure perturbation moving in upper layers of solar photosphere with supersonic speed is 
also considered. It is shown to excite acoustic waves putting certain restrictions upon their 
geometry of the generated wavefront. The results are discussed in relation to flare generated 
sunquakes. 
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1 INTRODUCTION 

Ray theory, as form of geometrical optics, has played a crucial 
role in helioseismology of acoustic oscillations, layin g foundations 



to su c h methods as time-distance helioseismo l ogy dPuvall et al 



19931: iKosovichev & Duval<l 19971: lGilej|200(]|: iD'Silva & Duval 



19951 : ID' SiM 19961: IXhompson & Zharkovll2008h. acous tic holog- 
raphy ('Lindsev & Braun"2000': lin dsev & BraunI |2004 and far- 
side imaging (Lindsev & Braun 200(J7 It has provided us with rich 
insight into the nature of acoustic wave propagation in the Sun from 
properties of low-degree modes and Duvall's law used in global he- 
lioseismology to averaging setups and focusing information relied 
upon in local helioseismic techniques. While more accurate numer- 
ical methods have since been deployed for computing helioseismic 
sensitivity kernels, which are used to deduce subsurface properties 
of the solar interior from observations, these have generally con- 
firmed and built on earlier results from ray-theory (ICouvidat et alj 
I2004l2006l) . 

Nonetheless, the use of ray-theory in solar physics has mostly 
centred on studies of individual rays, or modes, with little attention 
paid to fully transcribing the initial conditions responsible for gen- 
erating the family of rays that fully describes the wavefield. For ex- 
ample, 'Sogdan (1997) has compared the relationship between the 
modal and time-distance formulations of local helioseismology, but 
by own admission neglecting to take into account the initial value 
problem in terms of ray theory. One of the simplest of such cases 
is a spherical monochromatic point source, which is considered in 
this paper. 

When applied to standard wave-equation the mathematical 
method of geometric optics is known by other names such as De- 
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bye ray-theory, geometric acoustics, geometric seismics, WKB or 
semiclassical approximation leading to a fair amount of confusion. 
In solar physics some of the mo re rigorous applicati ons of the 
method can be found in works bv lOoughl Jl993L l200l and refer- 
ences therein). 

Here I use pu re ly mathematical approach as set out 
by iKravtsov & Orlovl ( Il99 0) for standard wave-equation and 
generalised for gener al asymptotic differential operators by 
iGuillemin & Sternberg! (tl990) to explicitly derive and interpret the 
eikonal and transport equations for non-linear wave equation in 
the presence of acoustic cut-off frequency. Finding the full solu- 
tion of the problem depends on reconstruction of the phase func- 
tion throughout the whole domain which is accomplished by solv- 
ing the eikonal equation. The phase function defines the geometry 
of the problem and can be found when proper initial conditions 
are prescribed. I then find such solution for the point source using 
basic polytrope model of the solar interior and consider its basic 
properties such as wavefront geometry, caustics and surface mani- 
festations. Q 

For physical applications, the results are compared with 
known acoustic sources in the Sun. Stochastic turbulent convec- 
tion in the upper layers of the convective zone is generally thought 
to generate most of the solar acoustic spectrum. However, another 
known source of sound waves in the Sun are acoustically active so- 
lar flares. In such flares, the energy release in the corona generates 
observable photospheric ripples that accelerate radially outward 
from a source region. This recently discovered phenomena known 



^ Note that there are also other independent means of deriving the 
eikonal equation such as Fourier transform and stationary phase method 



iChapman et alJl999l : lFelsen & Marcuvitj|l99ll) . 
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Figure 1. Spatial geometry of the problem 

as 'sunquake ' dKosovichev & Zharkovij|l998l) provides us with a 
localised example of wave excitation in the Sun. The acoustic na- 
ture of sunquake photospheric ripples has been well established and 
localisation of the source position can be deduced with good pre- 
cisio n using he lioseismic technique s such as time-distance anal- 
ysis dKosovichev & Zharkova 1998i:lKosovichev 2006. 2 7, ) an d 
acoustic holographylPonea et al.l ( lT999h : lDonea & LindsevlJioOsh : 
iLindsev & Doneal ( |200 8|). In particular, the egression computation 
in acoustic holography uses the Green's function built on the as- 
sumption of a monochromatic point source generating downward 
propagating waves from given depth in the solar atmosphere. 

The paper is organised as follows the mathematical formal- 
ism for asymptotical treatment and solution of the solar wave equa- 
tion, splitting it into eikonal and transport equations is described 
in Section [2] The link between general solution and initial condi- 
tions is addressed in|3] The solution for eikonal equation for a point 
source is derived in Section|4]including the expression for the Jaco- 
bian. Possible applications and limitations of the derived solutions 
to sun-quakes and basic properties of moving source wave-fronts 
are discussed in Section|5] 



2 EQUATION FOR SCALED PRESSURE 
PERTURBATION AND GEOMETRICAL 
ASYMPTOTICS 

Solar wave equation for scaled Lagrangian pressure perturbation, 
= /9q ^ Sp, can be written ( lGough|[l993h as 

Here Vh 's the horizontal Laplacian operator, po unperturbed quiet 
sun density, adiabatic sound speed, 7V^ Brunt- Vaasala buoyancy 
frequency, and 

„2 



C 



2 (1 - 2n. V H,) 



(2.2) 



is Lamb's acoustic cut-off frequency, defined in terms of density 
scale height Hp and unit vector n in the direction of the grav- 
ity action. As is standard in acoustic mode helioseismology let us 
take A^^ = to simplify the equation to a non-linear second-order 
wave-equation: 

4('£+^L')*-v'* = (2.3) 



C2 \dt^ 

To find asymptotic solutions the method of geometrical optics 



dKravtsov & Orlovll 19901) is used in this paper. The method, devel- 
oped for linear wave-e quation and generalised as met hod of Geo- 
metric Asymptotics bv iGuillemin & Sternberg! dl990l) for a wider 
range of high order differential operators. In short, this works by 
looking for the solutio n of the equ ation ( 12. 3t in terms of Debye ray 
series (or WKB ansatz IChapman et al.i tl99ft) ) written for asymp- 
totically large parameter A: 



(iA) 



(2.4) 



On substitution of the series into ( 12.3) and arranging the result in 
terms powers of the large asymptotic parameter A a series of equa- 
tions is obtained: 

2 2 
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c2 dt dt 
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(2.5a) 



(2.5b) 



2 dip dAi 
'c2 dt dt 



dt^ 



V 



c 



1 d^ 

■2 df2 



-— Uo (2.5c) 



Here ( I2.5at is a first order differential equation in partial 
derivatives called an eikonal equation which under suitable initial 
conditions can be solved to fully reconstruct phase function ip{r, t), 
which is also known as eikonal. The following J2.5bl2.5ct are 
called transport equations and can be solved iteratively using the 
solution for eikonal to obtain amplitude coefficients Ao, Ai, . . .. 
For more details and exam ples of the application of the method in 
the solar case see lGoughidl99l . 

Let A/ be a space-time manifold, described by Cartesian co- 
ordinates (x, y, z) and time, t. Let z be the depth, and x and y 
correspond to horizontal directions. Then, as <^ = ipix, y, z, t), the 
following variables describe ip in the phase space, Q,Ad 



_dip _dip _dip 

dx dy dz 



The eikonal ( |2.5a| > belongs to the class of Hamilton-Jacobi equa- 
tions and can be solved by method of characteristics in the phase 
space. In the above notation it corresponds to the Hamiltonian 



(2.7) 



so that the equation l |2.5a| ( becomes H — and finding the solution 
is reduced to integration of the characteristic system. 

It can be shown (see Section [X2l for example) that when the 
Hamiltonian ( |2.7| l does not explicitly depend on horizontal vari- 
ables, under suitably symmetric initial conditions, AI can be re- 
duced by one dimension, so that the solution can primarily be 
sought in two spatial dimensions, i.e. (p — ip>(x, z, t), so that 



1 f , 2 ,2 UJ ~ u 
H =-ikh + k, -j- 



dp , dip dip 
dx dz dt 



(2.8a) 



(2.8b) 



The full solution is then obtained via rotation using kf^ — k^ + ky 
relationship. 

In the problems considered in this paper, the sound speed 
and acoustic cut-off frequency depend only on depth, so the time- 
dependent equation i2.3\ is often reduced to two spatial dimensions 
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plus time. In such case we introduce cartesian coordinates x and 
z on manifold M, corresponding to horizontal and vertical direc- 
tions, with depth z chosen to be at the surface and positive below 
it (see Figure[T}. Note that in six dimensional space Q.M the graph 
of (i.e. the map (x, 2, t) M> (x, |f, is a three- 

dimensional surface. The characteristic system is then defined as 



dx _ dH_ dz _ dH_ dt _ dH 
dr dkh ' dr dkz ' dr dio ' 
dkh _ dH dk^ _ dH dm _ dH 
dr dx ' dr dz ^ dr dt ' 



(2.9a) 
(2.9b) 



where t is an independent variable. Note that condition H — 
can be used instead of one of the equations in the above sys- 
tem. Solutions of the system ( 12.91 ) written in terms of r (i.e. 
x(t), 2;(r), t(r), fch(r), ^^(r), cj(r)) provide us with characteris- 
tic lines lying on the graph of (f in Q,M. The phase function is 
determined by integrating 



dip dx dz dt 
1- = ^h-r + i^^'-i — ^-r 

dr dr dr dr 



(2.10) 



along each such line. In order to reconstruct the phase function 
if in the whole (a;, z, t) domain, initial conditions on some two- 
dimensional surface S G Q,M have to be provided so that charac- 
teristic lines originating from every point of such surface sweep out 
the three-dimensional graph. Let {C,, rf) be the coordinates on S so 
that the initial field can be written as 



(lA)'" 



Characteristic solutions of ( 12. 9t with initial conditions taken as de- 
scribed sweep out a three-dimensional surface in Q.M described 
by ((^, r), t) which are called ray-coordinates. Thus, it is impor- 
tant to take into account the proper initial conditions in order 
to understand the geometrical aspects of the wave-field associ- 
ated with the equation ( 12. Sat . Projections of individual characteris- 
tics (a;(r), z{t) ,t{T) ,kh{r) , kz{T),u}{T)) from Q.M onto {x,z,t) 
space are called rays. 

Once the eikonal equation solved, its solution can then be used 
to reconstruct amplitude along the rays by rewriting the transport 
equations in ray coordinates. These take form: 



dAl ^dWDij) 
^ + ^» dr 
^_dA, ^ ,^ d\nV{T) 
dr dr 



0, 



. , 1 a^Ao 



(2.11a) 
(2.11b) 



Here T){t) 



a(x,z.t) 
e(C,»7,T) 



is the Jacobian of the transform from 
cartesian to ray coordinates. Zero-order equation has the solution 
dKravtsov & Orlov.,199Q) : 



Ao{C,v,r)^A°{C,v) 



V{r) 



(2.12) 



Note that Aq depends only on geometry of the problem and the 
initial values of the field along the ray, thus expressing the local 
character of zero-order WKB approximation. So for example, in 
order to reconstruct the field at some point {x\, z\,t\), one nor- 
mally determines rays coordinates corresponding to that point via 
ekional equation solution and then uses ( 12. 12t to find the value of 
Aq and ( 12. lot to find the phase-function. 



Higher order coefficients can be written as 



V{r) 



V Ar, 



_ 1 d''Am-l 
' C2 dt^ 



^/V(T)dr'. (2.13) 



The presence of the Laplacian and higher order time derivative in 
the integral on the right hand side in the above means that Am de- 
pends not only on the values of Am~i immediately on the ray, but 
also the values Am-i in ray vicinity. Thus higher order amplitudes 
describe diffraction effects discarded in zero-order approximation. 

At points where the Jacobian is zero, Ao and higher or- 
der coefficients grow infinitely large, thus defining regions where 
the method breaks down. Such regions requir e separate analysis 
such as matching asymptotics or other methods iKravtsov & Orlovl 
1990; Kravtsov & Orlov 1993). Points where the Jacobian is zero 
are said to form caustic surfaces, or caustics. Such surfaces enve- 
lope ray surfaces and are singularities of projection o f the graph of 
fiC, rj, t) from QM to {x, z, t) space ( lAmolj|l978l) . 

It is common in many applications to disregard the non- 
zero order terms in ( 12.4) , Am = 0,Vm > 0, thus assuming 
the wave-equation ( I2.3t solution to be an "eikonal approxima- 
tion" 'i! Ki Aq exp(iAi/p) , also known a s an "almost plane wave". 
IKravtsov & Orlovl ( ll990l) : lGougri ( I2OO7I) give necessary conditions 
for this as that the medium, (? and oj^^ varies slowly over the char- 
acteristic length-scale of variation of the wave. 

3 INITIAL CONDITIONS 
3.1 General approach 

In 3 + 1 dimensions consider initial surface S C M parametrised 
by some coordinate set [6, (j>, to), with M described by {x, y, z, t). 
Then the field on such surface can be written as 



Solution of the 3D characteristic system 
{ dx dH 



dr 
dk 



dH 



(3.1) 

(3.2a) 
(3.2b) 



deduced from the Hamiltonian ( 12. 7t . provides the graph of (p, 

dip dip dp dp^ 



^5 ^! ^? 



dx' dy' dz' dt 



C Q.M, 



as function of ray-coordinates (9, (j>, to, r). Here r is the parameter 
along each ray, r = on S*, and 6, (j> and to are ray-numbering 
coordinates. As S is submerged in M, there is a natural mapping, 
{9, <f), to) ^ {x, y, z, t), providing the initial conditions for ( |3.2at . 
To recover the initial conditions for momentum part (13.2a) of the 
characteristic system, recall the definition ( 12.61 of the momentum 
variables, and let kxo, kyo, kzo,i^o be their values on the initial sur- 
face S. It then follows that 

' dpo dx dy dz dt 

^=fc^o^+fc.o^ + fc.o^-a;o^, (3.3a) 

dpo dx dy dz dt 

= ^'d^ ^ '"d^ ^"90 - "°a0' ^'-'"^ 

d^o dx dy dz dt 
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which together with 3-D Hamiltonian J2.7t evaluated at point 
{6, (ji,to, kxQ, kyo, kzOjOJo) forms a set of four equations for the 
four unknowns (kxo, kyo, kzo, oJo). These are then solved in terms 
of ray-numbering coordinates {6,4>,ta), thus recovering full ini- 
tial conditions for the characteristic system. Solution of the system 
( 13. 2t provides the map 



(61,0, to, t) h> (x,y,z,t,kx 



h h 



sweeping out the graph of the phase function 'p{x, y, z, t) in the 
phase space from the initial surface S. The phase function itself and 
amplitude coefficients are then reconstructed as described above. 



3.2 Monochromatic spherical point source 

Let us now derive the initial conditions for a spherical source. Con- 
sider a static sphere of a small radius r located at depth z = Zs 
to be the surface S where the initial conditions for pressure per- 
turbation are set. The surface can be parametrised by the spherical 
coordinates plus time, (6, cj>, to), so that the Cartesian coordinates 
and time are expressed as 



X = Xs + r cos sm <; 
y = ys + r cos 6 cos < 
z = Zs + r sin 6 
t = to, 



where (f) G [0, 2tt),9 e 



TV TT 
"2' li 



. Using ( I3.U3.3I >. let kl 



(3.4a) 
(3.4b) 
(3.4c) 
(3.4d) 

" k'io + 



kyo and define fc^o so that fc^^o = fcho cos fc^o, kyo = k^o sin k^o, 
system l l3.3t then simplifies to 



de 

difo 
d(j) 
dipo 
dto 



kzor cos 6 — khof sin 9 sin (fc^o + </>) , (3.5a) 
fcho'" cos 6^ cos (fc^o + <?!') , (3.5b) 



(3.5c) 



When sound speed and acoustic cut-off frequency depend only on 
depth, Hamiltonian J2.7I ) is explicitly independent of x, y and t, so 
from ( I3.2al3.2b] l it follows that kx , ky and uj are constant along 
each ray. Therefore, kh and k^ extended from the above definition 
throughout the phase space are also constant along rays. When, in 
addition, the initial phase function ipo does not depend on 0, e.g. 
9^50 /d4) = 0, it follows that cos (fc^o + 0) = 0, so that via a 
suitable choice of the x- and y-axis directions fc^o + <t> = + 2im. 
In this case the system can be reduced to two-dimensions plus time 
by letting = ^, fc^o = in l l3.4t . introducing x = x — Xs, and 
solving system J2.9b . i.e. using {6, to) to find (i, 2, t, fc^, fc^, cj). 
Using symmetry, this solution is then rotated by varying (j>, k4, to 
obtain the full three-dimensional field: 



(6,to,T) {x,z,t,kh,kz,Ld) ^ {x,y,z,t,kx 



(3.6a) 
(3.6b) 



X = Xa + X sin (l>,y = ys + X COS (p, 

kx = kh sin cj>, ky = kh cos (j>. 

Point source is approximated by letting r in the definition 0.4t 
to be infinitely small. A source is spatially homogenous when the 
initial phase function, ipo in ( 13.1b . may depend only on time. 

When the dimensionality is reduced, surface S, where initial 
conditions for the pressure perturbation are set, becomes a circle of 



radius, r, located at depth z 
coordinates {d, to) with 



The surface is parametrised by 



X = Xs + r cos o 
z = Zs + r sin 6 
t = to 

The initial field on this surface is described as 

■fo = Ao{e,to)e"^<'^''-'°\ 



(3.7a) 
(3.7b) 
(3.7c) 

(3.8) 



The unknown initial wavenumbers and frequency, kho = 

d(po/dx,kzo — difo/dz and ojo = —dipo/dt, are found from 
the system: 



dipo dx dz 

ato oto dto 

dipo _ , dx dz 

-Qf-'"^OQ^+kzOQ^- 



- dJO 



dt_ 

dto' 



dt 



khO + kzo — 



2 



(3.9a) 
(3.9b) 
(3.9c) 



Let us now consider spherical monochromatic homogeneous 
source of some fixed frequency, uuf, i.e. ipoid, to) = —ujjto. Then 
from the above ujo = ujf and kho sin 6 — kzo cos 6. This implies 
that in this configuration (see Figure[T) 9 can be viewed as ray take- 
off angle and all rays generated from S are of the same frequency, 
so frequency subscripts can be dropped. Define 



ks{Zs,U}) = 



then horizontal and vertical wavenumbers can be rewritten as 
kho = ks cos 9, and kzo ~ fcs sin 6'. Since Hamiltonian is inde- 
pendent of horizontal coordinate and time, lj and kh are constant 
on each ray. Then so is the horizontal phasespeed of a ray, 



2 

Vph 



'^h 



When to is fixed, 6 serves as a ray-numbering parameter. Note 
that by setting A{9, to) = 0, when 6 < one can investigate the 
wavefield generated by rays going down from the source, which 
would correspond to semi-spheric pump creating pressure at the 
source. Similarly, taking A{9, to) — for 6 > corresponds to 
the wavefield generated by rays going up to the surface. 



4 POLYTROPE 

4.1 Individual ray solution 

Let us now consider the reduced system in 2 + 1 dimensions with 
Hamiltonian J2.8at independent of horizontal coordinate and time. 
The characteristic equations ( 12. 9b are then considerably simplified: 



dx 
d7 



dz 
doj 



— kz 



dt _ uj kh _ Q — _ Q_ 
dr c^' dr ' dr ' 



kl + kl 



= 0. 



(4.1a) 
(4.1b) 
(4.1c) 



When kh = 0, the corresponding ray propagates purely vertically 
(x = const). It is easy to see that in terms of the reduction de- 
scribed for spherical source in Section |3!2l this ray travels along the 
rotation axis. 
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Consider polytrope model with adiabatic sound-speed and 
acoustic-cutoff frequency depending on deptfi only: 



C (z) = , u^^{z) 

m 



g(m + 2) 

42 ' 



(4.2) 



2.67 X 



where g is the gravitational constant (the value of ; 
10"'' Mm/s^ is used in our calculations), and m is the polytrope 
index. Let us discuss some general properties of such model and 
eikonal solving associated rays. 

A single ray with horizontal wavenumber kh and frequency 
iu in this model, in general, will have two turning points, Zu and 
zi (upper and lower), determined by the condition kz — 0. In the 
exceptional case of ray propagating purely vertically, kh = 0, the 
lower turning point can be viewed as located at infinite depth. Oth- 
erwise, using i4.2\ one obtains: 



A-2 



2 I 
'Z + 



m(m + 2) 



kl9 

Zu)(zi - 



-af 



(4.3) 



where Zu,zi are the upper and lower turning points and a[kh,oj) = 
{zu + zi)/2, b(kh, Lo) ~ (zi — Zu)/2. These can also be written as 



a{kh,uj) = 



1 cj m 
2^ 



,b(kh,uj) = , a'^{kh,uj) 



m{m + 2) 



Note that the upper and lower turning points coincide when 6 = 0. 

Furthermore, in this model the wave-vector length, k^ — 
(cj^ — uj'^^)/c^, viewed as function of z only, has a maximum at 
g{m + 2)/2a;^ determined by lj^^{ze) = uJ^ /2. This value 



ze 



is important due to the following more general statement: 



Proposition 1. With Hamiltonian H in ( I2.8at independent of x 
and t, every ray solving ( 12.9) that has two distinct turning point 
depths, Zu and zi, will pass through at least one point with z = ze, 
where ze is such that 



d_ 

dz 



- 0. 



Indeed, according to ( 12. 9b . 



dkz 
d7 



dH 

dz 



d_ 

dz 



. But this is the rate of change of vertical wavenumber along the 
ray, and since fcz(r) is smooth and becomes zero at turning points, 
it will have an extremum between them. 

Corollary 1. If for a given frequency there exists a unique "parti- 
tion depth" z = Ze such that 



dH 

dz 



0, 



then for all rays with distinct upper and lower turning points the 
following inequalities are true: Zu < ze,zi > ze - Ray starting at 
depth Ze with 



;,2 _ ^ 

— 



will propagate horizontally. 

In the polytrope model, b becomes zero only when ray is ini- 
tialised at depth ze{uj) with kh ~ ui^m/2gzE- It follows that un- 
der these conditions a = ze- 



Since cj/c^ 7^ in the domain of interest, z > 0, one can 



rewrite system ( 14.1) using group travel time, tg, as the parameter 
along the ray, with t — tg + to, to obtain 



dx 
dig 
dcj 
din 



kl 



kh, 



dz 

dta 



0,^ = 0, 
dtg 

■ b^ — {z — a)^ 



(4.4a) 
(4.4b) 
(4.4c) 



Therefore when 6 7^ 0, 



hence. 



2; = a + fe sin a 

1 CJ 



khg 



tg + Cz, 



X = 



2 kh 



-to 



b cos a + Cx, 



(4.5a) 
(4.5b) 



where a = tgkhg /ojm + Cz and integration constants Cz and 
Cx are determined from the initial conditions. For & = it is 
easily checked that z = ze = a,x = u)tg/2kh + Cx, i.e. 
( |4.5b| l still holds. In the case when kh = Q, x = const, z — 
gt^g/ 4m ±Cztg \/ g/m + Cl, where the choice of sign depends on 
the direction of the ray propagation and Cz is detemined from the 
initial conditions. 

Using group travel time as a parameter along the ray leads to a 
following interpretation. Let to) be the coordinates on the initial 
surface with the field described there as 

*o = ^E!(e,io)e'^'''«'*''', 

such that = when to < 0. Then from ( 12.12) it follows that 
Ao{(^,t(i,tg) = along the ray when t = to + tg < tg. Hence, 
group travel time is the time it takes for the initial perturbation to 
travel to a point along the ray and surfa ces tg = const represen t 
wavefront snapshots at a particular time dKravtsov & Orlovll 19901) . 

Using the monochromatic spherical source initial conditions 
( 13.713.8) with a;^ = the ray system can be solved in terms of co- 
ordinates {0, to,tg) as shown in Appendix lAl Given source depth, 
Zs, the solution is written in terms of take-off angle, 6, and group 
travel time tg- The group travel time is expressed via a parameter 
a{8, tg) and its initial value ao(0, Zs)- Solution for ao is given in 
( IA4) . which reveals the importance of the source position relative 
to the turning point partition depth, ze- Formulas dAlalAlet to- 
gether with the system ( |A5) provide the immersion of the phase 
function graph into phase space: {9,to,tg) ^IM- Formulas 
dAlalAldt themselves represent the projection of such map to M: 
(6, to, tg) n- {x, z, t)- Due to the simple translational relationship 
dAlc) between the group travel time and time and initial time vari- 
ables, when considering eikonal solution one is justified in simply 
fixing to- Then for 6 = const, JAlalAld) provide parametric ray 
solution as a function of tg, while taking tg = const gives geomet- 
ric wavefront at the time, t = tg + to, as function of take-off angle 

e. 

Examples of the solutions in terms of ray and wavefront prop- 
agation for monochromatic spherical sources placed at two differ- 
ent depths are presented in Figure |2] In both cases, the polytrope 
index is set as |,to — 0, and source frequency, 1/ = oj/2-k, is 
6.5 mHz. Each panel represents a spatial snapshot of the eikonal 
system at a particular time which is given in the panel title. The 
wavefront is plotted as thick red curve, selected ray paths up to the 
time indicated are plotted in black. The left column corresponds 
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Figure 2. Acoustic wavefront at different times generated by monocliTomatic point source with frequency u = 6.5 mHz with Zb = 1 Mm (left) and Zs = -142 
Mm (right). Rays start at to = 0, slice times are given in plot titles, geometric wavefront is plotted in red, with selected time-limited rays in black. Caustic 
points are marked by blue squares. The polytrope index is taken as m = | . The surface is located at z = 0. Note that the right column corresponds to 
the case when the source is placed in the region where variability of the model is high, so that the eikonal approximation becomes invalid. Indeed there are 
reasons to think that such a solution does not represent a physically realistic scenario (C. Lindsey - private communication). See Section IB. l.ll for discussion 
of applicability. 
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Figure 3. Downward propagating rays with u = 6.5 mHz, m = ^ for 
point-sources at deptlis of 1 Mm (blue) and 2 Mm (red). Envelopes of the 
ray families representing the caustics are clearly seen. 



to the case Zs — 1 Mm > ze, right column shows snapshots for 
Zs = 0.142 Mm < Ze- After simultaneously leaving the source, 
rays propagate up (9 < 0) and down (6 > 0) from the source 
sweeping out the generated wavefront. Near the surface, at the up- 
per turning points rays reflect back to the interior due to acoustic- 
cutoff frequency increase. In the interior the rays are reflected back 
to the surface at lower turning points. As different rays reach turn- 
ing points at different times, the wavefront gets deformed. This is 
discussed in the next Section. When source is located below the par- 
tition depth Ze (left column), wavefront reaches the surface with 
the rays going up from the source and then propagates radially 
outward in the near-surface layers with the second bounce surface 
manifestation appearing later. When source is located above ze, 
the situation is different: after initial surface bounce for rays going 
up from the source, the reflected wavefront goes below the surface 
mirroring the behaviour of the rays going down, before reappear- 
ing some distance away from the source and then travelling radially 
outward in the near surface layers. This is further discussed in Sec- 
tion|43] 



4.2 Caustics 

Figure |3] shows two families of rays generated by semi-spherical 
point sources located at 1 and 2 Mm depths. In both cases, the 
downward pump scenaricQis represented, i.e. Aq{9, to) — when 
{6 ^ 0), so only downward propagating rays are plotted. It can 
be seen that near the upper and the consecutive lower turning 
points the density of the projected onto M characteristic curves (i.e. 
rays) increases with ray families becoming enveloped by "caus- 
tics". Caustics are defined as suifaces where the Jacobian of the 
coordinate transform from ray to Cartesian coordinates becomes 
zero. This implies that at such points the mapping is not unique, 
i.e. two or more rays pass through the same point in {x, z, t) space. 
It can be checked that, under the initial conditions considered here, 
this is equivalent to the wave-vector being tangent to the wave- 
front. As wavefront propagates in the direction of it's tangent and 
rays travel in the direction of the wave-vector, it follows that such 
points locally envelop the wavefront and ray-paths. Caustics are 
often, though not always, associated with the singularities of the 



^ see the last paragraph of Section l3.2l 



generated wavefronts su ch as cusps. See dKravtsov & Orlovll 19901 : 
iKravtsov & Orlovl 19931) and references therein for more details. 

Fore G (--J, -j) , analytical expression of the Jacobian of the 
transform from ray-to-Cartesian coordinates is derived in Appendix 
I A2I (see JAlOb ). Evaluating J (6, a) at lower turning points, a = 
■| -I- 2n7r, gives J{9, ^ + 2mv) — —zftand/b, and at the upper: 
J{6, -f + 2n7r) = zltanO/b. It follows then that the turning 
points are on caustics only when ^ = 0. Moreover, 

J(0, a) = bzs cos a ^1 - = -2b^-^cosa, 

so that the ray with 6* = touches caustics only at its turning 
points. In particular, this means that the initial value of the Jaco- 
bian is zero, which according to ( |2.I2t indicates that in zero-order 
approximation the amplitude of the wave-filed will be zero along 
the ray. When take-off angle is non-zero, the fact that the Jacobian 
along the ray takes different signs at upper and lower turning points 
suggests that there is at least one caustic point situated in between. 

The caustic points, computed numerically using dAlOt . at 
given times are plotted as blue squares in Figure |2] In both cases 
caustics are formed after rays went through the turning points. 
When Zs < ze, right column, the first caustic is formed at the 
interface of the downward propagating front and the reflected from 
the surface following front. Once this caustic comes back to the 
suiface (the bottom plot), it then propagates near the surface away 
from the source. In the case when Zs > ze, left column, the first 
caustic point is formed near the surface when the ray with 9 = 
passes through the upper turning point. The point then travels ra- 
dially in the near-surface layers away from the source. The bottom 
image shows the appearance of the next set of caustic points, one 
of which is again located near the surface and is associated with 
the second bounce waves at the surface. The others correspond to 
wavefront cusp formed near lower turning points. A close-up ex- 
ample of the formation of caustics associated with the wave-front 
cusp is shown in the Figure|4] 

In Figure |5] locations of the caustics and turning points are 
compared. It can be seen that the lower caustics are associated with 
lower turning points for second and further bounces. These caustics 
are quite close to the locations of the turning point, though not ex- 
act except for 6^ = 0. The caustics formed near the surface are very 
close to rays' upper turning points. Moreover, near surface caus- 
tics are generally quite basic, e.g. not associated with wavefront 
singularities, and are formed when only two rays (one after being 
reflected and another coming up to its reflection point) intersect. 

The results presented in this Section agree well with the anal- 
ysis of the numer ical simulation of a coustic wave propagation in 
the solar interior dShelvag et al.ll2009l) . Further analysis is needed 
to rigorously determine types of the caustics generated by a spheri- 
cal point source and hence describe acoustic field behaviour in their 
neighbourhood. 

4.3 Wavefleld at the surface 

As no direct observations of the Sun below the photosphere are 
available, the surface manifestation would be a primary observ- 
able for any acoustic wave-field generated in the interior or near 
the surface. Let Xsurf{9\ijJ, Zs) = Xsurf{9) denote the horizontal 
distance from the source to the first appearance of the ray on the 
surface, i.e. its first upper turning point. As shown in Appendix lAll 
upper turning points correspond to a = + 27m, so the first 
such point away from the source depends on the choice of initial 
value, ao- 
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Figure 4. Formation of the wavefront singularity and caustics near lower turning points, for Zs = 1 Mm, = 6.5 mHz. The caustic is formed when "up" rays 
after being reflected near surface, intersect "down" rays travelling up to the surface after passing the lower turning point. 



When source is placed at or below ze, Zs ^ ze, the situation 
is straight forward with Qo(6l) G (-f , ^) for 6^ G (-f, f) , 
so that the first upper turning point along every ray corresponds to 



^TT. Then 



surf (^) — 



where 



■2cos2 6l V2 



TT — Qfo + 2s tan t 



(4.6) 



(see Appendix lAll for more details). The corresponding group 
travel time, tg,surf, the time it takes for the perturbation to travel 
along the ray from the source to ray's first upper turning point, is 
given as 



g .surf 



khg 



- ao 



(4.7) 



When source is placed above ze, i.e. Zs < ze, ao is less or equal to 
— ^ when 9 is negative (rays going up from the source), and greater 
than or equal to — .| when 6 is positive. Thus, 



^ sur f 



'1 _e 

^ 2 cos^ 

'1 « 



k 2 cos-i 



) + tan 9) 
+ tan e) 



9 <0; 
9>0, 



with corresponding group travel times: 



tg,surf{9) — 



^(-i^-ao),e<0; 
^(|7r-ao),e^O. 



Due to the symmetry of the initial conditions, the surface rip- 
ples will propagate circularly away from the source, following the 
time-distance relationship defined above using the take-off angle 
as parameter. Comparing these time-distance relations with the one 
for skip-distance, A = t\g/^-Km (see Appendix lAll for details), 
it is useful to separately consider rays going up and down from the 
source. The rays going up reach their upper turning point first and 



reflect back to the surface. After that the first surface ripple is de- 
fined by the rays leaving the source away from the surface and then 
reflected back. 

• when Zs ze, after the up-propagating rays are reflected 
from the surface, the time-distance curve is similar to the skip- 
distance relationship with times slightly scaled, since the rays 
do not travel the whole skip-distance (distance between surface 
bounces). This scaling/difference depends on the source frequency 
and depth, via parameter ^. See left column of Figure|2]for an ex- 
ample of such source. 

• when Zs < ze, after the initial surface bounce of "up" (9 < 
0) rays, the time-distance curve is even closer to the skip-distance 
relationship, with wavefront from down rays slightly pre-ceding the 
skip-distance arrival time, while the up-ray wavefront which was 
initially reflected from the surface near the source, arriving slightly 
later due to the time it takes for the ray to travel from the source to 
the surface. See right column of Figure |2]for an example of such 



5 DISCUSSION AND APPLICATIONS TO THE SUN 

Before moving to the discussion of the applications of this method 
to the Sun, let us briefly review existing observations of the acoustic 
wave-fronts. While acoustic oscillations in the Sun have been ex- 
tensively studied via global and local helioseismology, our under- 
standing of the nature of their excitation remains relatively sketchy 
with stochastic turbulent convection generally considered to be the 
main mechanism responsible. There is, however, a recently dis- 
covered sunquake phenomena that provides us with observational 
examples of locally generated acous tic fields. Sun quakes, which 
are relatively rare events (please see iDoneallioT l[ for discussion 
on detectability), take place in active regions during flares, when 
the flare induced changes in sunspot penumbra under the right cir- 
cumstances generate acoustic waves. These waves can be seen in 
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Figure 5. Caustics and turning points: projection of caustics onto spatial 
plane are thick black lines. Lower turning points of rays going down from 
the source are drawn in red, with lower turning points for rays going up 
from the source drawn in green. Upper turning points are in purple. 



photosp heric Dopplersram observations as circular shaped surface 
ripples jKosovichevI |2004 iMoradi et al.] l2007l : iKosoviche-vl 1201 ih 
accelerating away from the source region. However, due to back- 
ground oscillations, such ripples are often difficult to distinguish, 
so helioseismic methods such as time-distance diagram analysis 
and acoustic holography are usually applied to suitably processed 
series of line of sight velocity ob servations for quake detec ti on. 

Time-d i stance met h od teosovichev & Zharkoval 1 19981 : 
iKosovichevI |2007| , 1201 ll : [zharkov et al. 201^ provides direct 
observational evidence of generated acoustic waves by integrating 
the signal over circles or arcs centred at the source for every obser- 
vation in the time-series, thus producing a 2-D diagram displaying 
the field dependence on time and distance from the source. In 
such diagrams quake events have a clear signature in the form of a 
discernible time-distance ridge. It is the close agreement between 
the detected ridges and the theoretical travel-time skip-distance 
relation that provided the conclusive proof of sunquakes acoustic 
nature. On the othe r hand, the acoustic egres s ion measurements 
jPoneaetalJ 1 19991: loonea & LindsevI l2005l : iLindsev & Done j 
boost [Matthews et al.ll20 111) are based on theoretical modelling 
of acoustic waves propagating from a point source and provide 
a proxy measurement of acoustic energy emitted at a location 
at given time, thus producing a map of acoustic sources and 
sinks. Quake signatures in egression power maps are normally 
represented by a compact kernel or kernels of enhanced emission 
surround ed by the acoustically absorbing sunspot interior (see 
[zharkov et al.l ( l20lT b. 2012) for discussion and comparison of the 
two methods). 

Both methods were originally based on the assumption of 
photospheric manifestation of a flare generated acoustic field in 
the form of circular shaped ripples propagating away from the 
source. However, considerable anisotropy in the acoustic ampli- 
tude of the ripples from t he vantage of the s ources has been ob- 
served for m ost quakes jKosovichevI l200d : iMoradi et alj |2007| : 
lDoneall20T l|) with acoustic emission stronger in some directions 
than others and ripples changing shape becoming elliptical Q In 
egression power maps such anisotropy is represented by "stacked 



acoustic kernels", where tw o or three narrow kernels a re stacked 
together. IPonea et aO ( Il999l) : IPonea & LindsevI | |2005|) suggested 
that this could be the result of interference caused by the rapid mo- 
tion of the source, roughly in the direction along which the kernels 
are stacked. This was supported by observations of apparent rapid 
photospheric movements at the quake locations. iDonea & LindsevI 
(2005) found that for 2003 October 28 and 29 quakes the motion 
of the HXR sources was indeed aligned accordingly with egres- 
sion po wer stacks. This w as further confirmed for 2002 July 23 
flare bv IKosovichevI ^20071) using HXR and Doppler data, where 
via data analysis and modelling the author estimated the speed 
of the seismic s ource to be s uperso nic|3 around 20 — 25 km s~^. 
More recently, KosovichevI ( 1201 ih used time-distance diagram 
analysis to detect a supersonic movement of one of the seismic 
sources produced by 2011 February 15 flare, with speed around 
15 — 17 km s~^. Supersonic motions of around 14 — 22 km s^^ 
have also been detec ted for this flare's s econd seismic source in 
IZharkov et al. I( l2012l) . In the recent review IPoneal (1201 1) notes that 
the maximum amplitude of the ripples emanating from a moving 
source is generally along the axis of the source, displaced from the 
source location in the direction of the motion. 

Both helioseismic methods have shown very local spatial na- 
ture of the seismic sources validating the point source assumption. 
However, there is no current consensus in regard to the physi- 
cal mechanism and processes behind the quake excitation. Back- 
warming heating, hydrodynamic shocks, particle precipitation and 
Lorentz force are the main scenarios currently considered and de- 
bated asthose capable of producing flare acoustic response (see 
iDonej l lioilh and references therein for more details). Due to the 
lack of our current understanding of excitation, little is known about 
the depth of the source, but it is generally assumed to be near the 
photosphere since various modelling methods show only a small 
fraction of the energy initially invested in the shock penetrate into 
the photosphere. 



5.1 Stationary monochromatic source 

Back to the model considered in Sections l3l4l it is clear that non- 
evanescent acoustic waves are present in the model if and only if 
^ 0. Then the first condition on monochromatic source be- 
ing able to generate acoustic wavefield is uj^ > uj'aci^s)- Let 
Zcuti'^) denote the point such that Ldac{zcut) ~ oj'^- Furthermore, 
as shown in Section 1431 the presence of a unique extrema in the 
k^{z) = (oi^ — ujac)/c^ warrants the existence of the (upper 
and lower) turning point partition depth ze- When such a parti- 
tion depth exists, placing a source above it (but with uj > LJac(zs)) 
will create a gap in ripples seen at the surface. Indeed, using the 
fact that for a fixed frequency the skip-distance for a particular ray 
is a function of its horizontal wavespeed, one can evaluate 



= k (zs), 



^ for example, using directional time-distance diagrams iKosovichevj 
i201lh has reported strong travel-time anisotropy in different arcs of the 
wave-front generated by the strongest seismic source of 201 1 February 15 
flare 



therefore ^ jiuP" — cj^c)- As uj\^ grows near the surface, 
z approaches Zcut and the phase speed grows infinitely large. As 
upper turning depths of the rays are greater than ze, and since near 
upper turning points ray propagate nearly vertically, it is clear that 
rays going up from point source will not travel far from the source 
horizontally before being turned back down into the interior. For 

* in the quiet Sun adiabatic soundspeed at the photospheric level is esti- 
mated around 7 — 8 km s^^ in the quiet Sun 
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rays going down from the source the horizontal distance from the 
source to its first appearance at the surface will be slightly shorter 
than its skip-distance, so that the first appearance at the surface will 
be approximately at A(min v^^). Thus, for example, for polytrope 



model considered in Section HTTI taking m — ^, at z — 122 km 
only waves with frequency v — uj /2n ^ 7 mHz can be generated. 
Moreover, given Zs, if (jJ is close to the value of the ujac{zs), the 
distance where generated waves surface away from the source can 
grow very large due to the fact that ^ {lj^ — ijj^^)I<? . 

The assumption of a unique extrema for fc'^(z) doesn't h old 
in the case of Model C tehriste nsen-Dalsgaard et al] Il996l) . a 
widely used and thoroughly tested in helioseismology model of 
non-magnetic solar interior. This is due to multiple peaks in acous- 
tic cut-off frequency present near the surface (see bottom panel of 
Figure [6). This effect, however, is thought to be numeric in na- 
ture (see ISchunker & Callvll2006 i for example) and other models 
such as isothermal cut-off frequency, (? /AHp, are often used in 
the acoustic oscillation studies in order to ensure smooth variation. 
It can be checked that in the isothermal case such an extrema exists, 
and corresponds to a maximum that appears to be unique below the 
surface for frequencies of up to « 7 mHz. 

Adiabatic sound speed and acoustic cut-off frequency esti- 
mated from this model are presented in Figure |6] It is seen from 
Figure |6] where sound speed and acoustic cut-off frequency esti- 
mated from the model are plotted, that in the layers where ujac is 
large c is of the order of 7-8 kilometres per second. Note that in the 
Sun high frequency rays will not necessarily have an upper turning 
point, escaping into the outer layers of solar atmosphere. 



5.1.1 Applicability of the method 

Before going further, let us note some considerable limitations of 
the model and solution using the shallow source considered in the 
earlier sections as an example. Stationary monochromatic source 
solution suggests that when the source is placed close to the surface 
the generated wavefront can surface some distance away from the 
source (right column, Figure |2). If this were to be applicable to 
sunquakes, it would imply that the detection of such events may 
require a field of view larger than currently used. And even then the 
detectability of such waves via local helioseismic methods might be 
difficult due to reduced wave amplitude at large distances from the 
source. 

However, clearly the model used here is too basic to expect a 
realistic representation of the physical phenomena. Apart from the 
glaring lack of magnetic field which we know is present at sun- 
quake locations, the Brunt Vassala frequency is a lso zero, A'^'^ = 
in (12. It . Employing ray-mechanical formalism iBames & Callvl 
1 I2OOII) have shown the importance of Brunt Vassala term for waves 
near the acoustic cut-off boundary. 

There are also limitations associated with the solution itself. 
Firstly, by applying the condition fc^ ^ and thus rejecting com- 
plex valued wave-numbers , the phenomena such as jacket modes is 
omitted, whereas those are needed to form a complete set of wave- 
equation solutions ( iBogdan & Callvi[l995l) . This can possibly be 
addressed via extending the solution to complex rays, and for ex- 
ample considering the signal from downward propagating evanes- 
cent modes, i.e. with imaginary fcz, and evaluating if and how it can 
transform into pure acoustic rays at larger depths. 

Furthermore, the dispersion relation, corresponding Hamilto- 
nian and rays associated with initial conditions are derived from the 
eikonal equation l |2.5at . The actual wave-field is then reconstructed 
by solving the set of transport equations ( I2.5bll23c] ...) and postu- 



lating or establishing that the resulting series (12.4) are asymptotic to 
the exact solut ion. Under certain condit i ons such as sl ow-variability 
of the model ('Kravtsov & Orlovlll99(]| ; lGougril2007h , an "almost" 
plane-wave (or eikonal) approximation can be justified, where the 
resulting field is approximated via its zero-order amplitude. The 
usual interpretation of rays as paths along which the energy prop- 
agates is based on such an approximation. This is clearly not the 
case for the shallow source where, due to the properties of acous- 
tic cut-off frequency, model variability becomes high. As sufficient 
conditions for the "almost" plane wave approximation are clearly 
not met, this needs further investigation, e.g. solving higher order 
transport equations or approaching the problem via different meth- 
ods.. In fact, based on more realistic wave modelling, there are rea- 
sons to think that physical situation is different. (C.Lindsey - private 
communication) . 

Thus, the properties of the eikonal solution discussed here and 
in the following section should be considered as hypothetical, pos- 
sibly posing questions about real physical scenarios and applicabil- 
ity of the geometric asymptotics method. 

5.2 Moving source 

Let us now consider a moving source. This is achieved by going 
back to ( |3.4| l and letting the source coordinates {xs,ys,Zs) de- 
pend on time to. It is clear that the symmetry is broken unless the 
source moves along the axis of rotation, so the three-dimensional 
equations are considered. When solving the system ( 13.3) , use the 
stationary source case outlined in Section 13.21 As dtpo /dO and 
d(po/d(f> are actually the same, kho and fc^o are defined in the 
similar manner, and with equations ( |3.5a|3.5bt . it also follows that 
cos (fc0o + </>) = 0. It is only the frequency definition ( I3.5ct that is 
affected by source movement: 

dtpo , , dxs 

— — = ko.v - a;o, where v = — — , 
ato dto 

represents source velocity. Similar to Section 13.21 let the initial 
phase in the definition l l3.8t be monochromatic and homogenous, 
i.e. ifo = —ujfto. Then, the initial ray frequency is given as 



LOO 



■ u]f + ko.v. 



(5.1) 



Hence, in this case, the source frequency is modified depending on 
the wave-vector: ray frequency which is constant along each ray 
varies from one ray to another according \o ujo = uj f ~[- kov cos 7, 
where 7 is the angle between the wave vector and source velocity. 
The frequency increases along the vector of the source movement, 
and decreases when the wave-vector points in the opposite direc- 
tion. 

Similar to stationary source (see Section |3!2i , it follows from 
the first two equations in ( 13. 3t that k^ ~ ko cos 6, kzo = ko sin 9, 
where fco 



(^0 



„)/c . Then use the dispersion relation H • 



0, where Hamiltonian H is defined in i2.ll . to find ko. Substituting 
( 15. U into the last expression, the following quadratic equation with 
respect to fco is obtained: 



kl{u'^ - 1) + 2^uko + 
c 



= 0, 



(5.2) 



where u = vcos{j/c). Consider the case when, as before, the 
source with frequency greater than the acoustic cut-off value, oj/ ^ 



LOac{zs{to)) moves with subsonic speed, v < c 
solution is 



< 1. The 



1 u ± \/!d 

ko± = - 5 — 
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where D = cj^ — uj^,. + uj^i.u . Clearly, under our assumptions 
D > Q. Since fco is defined as length of the wave-vector, we are 
looking for positively valued roots only, which leaves only fco+ and 

To visualise the nature of the angle 7 let us now define 
an appropriate coordinate frame. Assume that at the time to the 
source is located at Xs{to) = ysito) = 0. Since there is free- 
dom in selecting the direction of x- and y- axes, choose coordi- 
nate axes so that the source moves in the direction of y-axis, i.e. 
V = {vx, 0, Vz). In spherical coordinates, (r, 6, (j>), defined for the 
source, v = {v, A,0), where A is the angle between the velocity 
vector and xy-plane. Then, using the spherical law of cosines, for 
any wave- vector ko = {ko,6,(j>) we have 



cos 7 = sin 6 sin A + cos 6 cos A cos < 



(5.3) 



Note that apart from the assumption of dependence on depth only 
all of the above does not depend on the model for the sound speed 
and acoustic cut-off. 

In the polytrope model the ray solution derived in the Ap- 
pendix, dAlalAlet together with ( IA4t and ( IA5t . holds for indi- 
vidual rays. However, the generated wave-field will have a more 
complex dependence on the parameter to, due to the source move- 
ment. E.g. at moment to source produces a family rays governed by 
the above equations, then at moment to + 5t the source will have 
moved producing another family of rays. This will be investigated 
in future publication. 

5.3 Supersonic source 

While the wave-fields generated by a monochromatic point source 
have some of the properties associated with solar quakes, one of the 
main limitations of the considerations above is the assumption of 
source frequency. Even for observed high frequencjjf] waves gener- 
ated by flares, the period is still of the order of couple of minutes. 
At the same time, the photospheric changes observed at t he quake 
locations happen on a much shorter scale feh arkova & Zharkovl 
l2007l : lKosoviche^^l200l[20Tll : IZharkov et alBoi la . l201 

It is instructive to consider a non-harmonic source, i.e. 
ifioiO, (j>, to) = 0. This corresponds to the wave-field on initial sur- 
face simply written as \E'o = Ao{6, (p, to), meaning that the pres- 
sure perturbation, p 2 Sp, is varying slowly on the initial surface. As 
in previous Section 15.21 dipo /dO — 0, fczo — kho tan 6, system 
i3.3\ leads to the following solution: 



ojo = ko.v = fco«cos7, 
,fcho ~ fco cos fczo = fco sin 6^, 



(5.4a) 

(5.4b) 
(5.4c) 



where ko is the wavevector, v is source velocity, v — |v|, and 7 is 
the angle between the two vectors. From the above it follows that 
non-evanescent acoustic waves (fco > 0) are generated if and only 
if the source moves with supersonic speed > c^. Moreover, let 
us rewrite ( I5.4bt as 



cos 7 



(5.5) 



where the square root is always taken with plus sign due to ( |5.4at 
and the fact that fco > by definition. Again, it is clear that while 

^ sunquake egression signal is usually strong around 1/ = 6 mHz 



the frequency is constant on each individual ray, it will vary from 
ray to ray as parameterised hy 6,(j> and to. It is also evident that in 
this case uo is always greater then the value of ujac at the source 



depth, going to infinity as cos 

^0 ^ ^min ~ 



7 



c /v < 1. In addition. 



- ujacV / {\/v'^ — c^). Note that when cos 7 is negative 
ujQ becomes negative, so waves waves will only be generated in the 
direction of the source movement. Moreover, given the (15.31 , the 
condition 



2 ^ c 
cos 7 > ^ 



(5.6) 



ensures that rays are generated for only a relatively narrow range 
of values of 6 and essentially forming a cone around the velocity 
vector V. The horizontal phase speed at the source does not depend 
on source depth: Vph ~ cjo/fcfeo ~ fcos7/cos^. 

Again, note that apart from the assumption of dependence on 
depth only all considerations in this section are independent of the 
model of the media where the waves propagate and rely only the 
dispersion relation as well as spatial geometry of the problem. Also, 
note that this solution is only possible when iUac > 0. 



5.4 Comparison with sunquakes 

As quake observations and modelling suggest the source is located 
near surface, let us assume that the source moves in the upper 
ranges of solar interior, so that for generated frequencies Zs < ze, 
i.e. the source is located near upper turning points. Then, for each 
generated ray, the horizontal distance from the source approxi- 
mately equals to the skip-distance, A{kh, uj). On the other hand, 
the skip-distance is essentially a function of horizontal phase speed 
of the ray with lower phase-speed values corresponding to smaller 
distances. The relationshi p is exact fo r polytrope and other theo- 
retical models (see Christe nsen-Dalsg aard 2003) and has been ob- 
servationally validated by the time-distance helioseismology and 
acoustic holography. 

Let us consider a source propagating with velocity v vertically 
downward, i.e. A = |- (see ( 15. 3t ). Then cos 7 = sin 9, and as only 
downward propagating rays are generated, the first appearance of 
the generated wavefront on the surface corresponds to the mini- 
mum value of Vph = wtanS. Then for 9 £ [O, ^] the minimum 
horizontal phase speed will be achieved at the lowest value of 9. 
Using l l5.6b one obtains 



mm Vph = 



(5.7) 



where expression on the right hand side is evaluated at z — Zs- 
However, frequency of the rays travelling near such phase speed 
will be approaching infinity due to ( 15. 5t . Hence if observations 
are made at certain Nyquist frequency, ojn, such waves may not 
be observed, so further restrictions need be considered, namely, 
ojo ^ tiJiv has to hold. From ( 15. 5t it follows that this condition 
is equivalent to 



sin 9 > 



From this the minimum "observable" phase-speed can be evalu- 
ated: 



obs 

mm Vph = 



(5.8) 
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Therefore, in the acoustic wave-field excited by a vertical super- 
sonic shock perturbation only waves with phase speed exceeding 
these will be observable. Moreover, as surface ripples from such 
source are determined by the phase speed, the minimum distance 
away from the source can be estimated using the above inequal- 
ity. For example, if c = 8km/s,i; = lOkm/s, then minimum 
Vph ~ 13.3km/s. For the Sun this corresponds to a skip-distance 
of around 5 — 7 Mm. Let us take the value of cut-off frequency 
at source depth, iUac{za)/'2n, to be ~ 5. mHz. Then adding the 
condition that the cyclic frequency along the ray is no greater than 
8.4 mHz, gives us a minimum phase speed estimate at ~ 106km/s. 
Therefore, in this case, only ripples at distances of order of hundred 
megameters from the source would be potentially observable. 
More generally, let us rewrite (15.3) : 

cos 7 = cos(^ — A) — cos cos A(l — cos 0). 

Then, as cos 9 and cos A are non-negative, the inequality cos 7 > 
c/v implies that cos{9 — A) > c/v > 0. Hence, 

^ G ( A — arccos — , A + arccos — ) . 
V V V J 

It is also clear cos > c/v — sin 9 sin A. Thus for downward prop- 
agating source, i.e. A G [O, -l] , very roughly it can be estimated 
cos (f) > c/v — sin A, so only waves in the limited range of (f) can be 
generated. But as rays propagate in the plane defined by z-axis and 
their initial wave-vector, this means that surface ripples will only 
appear in limited arcs in the direction of the source movement in 
this way producing anisotropy in wavefront amplitude. 

Given supersonic movements observed at the sunquake source 
locations, and strong anisotropy in amplitude in the direction of the 
movement, it is tempting to hypothesise that this mechanism plays 
an important role in the sun-quake wavefront generation. However, 
caution needs to be exercised as the problem has been approached 
here from purely mathematical standpoint and the notoriously diffi- 
cult question of method applicability has not been addressed. While 
this solution appears to provide a possible explanation for some of 
the observed properties of the sun-quake generated acoustic fields, 
the applicability of the method under the circumstances should be 
questioned as the variations are obviously very fast and known 
"sufficient" conditions for method applicability do not hold in this 
case (see Section [5.1.1| l. Numerical simulation or modal approach 
can provide a way to check this. It is interesting to note that pub- 
lished simulatio ns showing generation of acoustic waves by con- 
vective vortices ( Kitiash vili et alj|201ll : IMoII et alj|201lh also re- 
port supersonic movements at such vortices. Clearly further analy- 
sis is required to address these issues as well as to investigate vari- 
ous dependencies that in this paper were touched upon only briefly. 

The anisotropy of flare generated acoustic wave-fronts has so 
far been attributed to the effects of magnetic field and more general 
plasma subsurface properties in the Active Regions where quakes 
take place. It is indisputable that such conditions affect and likely 
shape generated waves. However, it can be argued that the effects 
discussed in this paper, if supported by more realistic analysis as 
well as for magnetised case, can play an important role in deter- 
mining the shape and properties of the wave-field, which is then 
inadvertently affected by the magnetic field, flows and other inho- 
mogeneities along its path. 

On the other hand, it can be seen that even without the mag- 
netic field, e.g. when moving source is considered, the number of 
variables on which solution depends grows quite fast. Given our un- 
ders tanding of interaction of acoustic and magneto -acoustic waves, 
(see ISchunker & cih] I2OO6I : iMoradi et al.l I2OI0I and references 



therein), it is clear that the inclusion of magnetic field effects is 
highly likely to paint an even more complex picture. 



6 CONCLUSIONS 

While the linear wave -equation has been extensively studied using 
geometric optics (see iKravtsov & Orlovlll990l ; iKravtsov & Orlovl 
1993, and references therein, for example), the non-linearity in- 
duced by the presence of acoustic-cutoff frequency in the Sun has 
often been overlooked. Using a point source as an example, in this 
paper I have shown how suitable description of initial conditions 
can be used to derive geometric properties of acoustic fields gov- 
erned by such non-linear wave-equation. For stationary monochro- 
matic source a full solution of the eikonal equation including Jaco- 
bian computation and phase function reconstruction was derived, 
showing dependency of the field shape on source depth and fre- 
quency. The general treatment of a moving source has also been 
outlined, analytically deriving some properties that might be as- 
sociated with sun-quakes. While primarily concentrating on math- 
ematical treatment to the problem, some intriguing properties of 
phase functions corresponding to near surface sources have been 
deduced. It was shown that due to the presence of acoustic cutoff 
frequency a shock travelling with supersonic speed can generate a 
cone-like wave-packet in the direction of the shock movement. 

Given the asymptotic nature of the method used, the paper 
poses interesting and current questions about acoustic waves gen- 
erated near the solar surface. The resolution of such questions will 
also shed further light on limits of applicability of the geometric 
asymptotics in solar applications. 

The availability of the full solution including the Jacobian 
suggests a possibility of using it to compute the Green's func- 
tion via solving of zero-order transport equation. This, how- 
ever, is complicated by the presence of the caustics, which 
require separate t reatm e nt and classificati o n as outlined in 
IKravtsov & Orlovl jl990l) : IKravtsov & Orlovl ( Il993l) . While this 
is a subject of future study, other methods su c h as outlined 
in 'P erez Hernandez & Gonzalez-HemandezI ( I2OI0I) : iLindsev et alj 
(2011, and references therein) might be better employed for this 
purpose. 

Nonetheless, the method presented in this paper is shown to be 
a powerful tool for studying the behaviour of the wavefields when 
full initial conditions can be specified. 
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APPENDIX A: MONOCHROMATIC POINT SOURCE: 
EIKONAL SOLUTION FOR POLYTROPE MODEL 

Al Characteristics 

In this section the eikonal equation for polytrope model is solved 
using the monochromatic point source initial conditions deduced 
in Section 1^21 For given frequency, uj, and source depth, Zg, we 
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reconstruct the graph of the phase function, ip, throughout the phase 
space domain in terms of ray-coordinates using group travel-time 
as the parameter along the ray. The Jacobian of the transform from 
ray to Cartesian coordinates is evaluated, and the phase function is 
reconstructed. 

Without loss of generality, take Xs = and cj > 0. By con- 



struction (see Section|3]2) , kh ^ and 6* G (-f , f ) ■ Note that 
by our definition (see Figure [T} 6 is zero in the horizontal direction 
and is increasing clock- wise. For fixed Za,Ljj, the initial values for 
wavevector components are determined in terms of the take-off an- 
gle, ^, and given by fcho ~ ks{zs,Ld) cos 0, k^o = ks{zs,uj)sm9. 
Using the corresponding upper and lower turning points, in accor- 
dance with ( 14.3b . define 



a{6\oj, Zs 



2 kig 



(1 + tan 



and 



,2,^1 \ 2 m(m + 2) 
b {e\uj,zs) = a . 



Then, in view of l l4.5bt . the solution can be written as 



la;,, , 

E = tn — (cos a — cos an) 

2kh 

z — a + fesina, 
t = tg + to, 
kug. 



a{tg,0\u], Zs) 
ao{0\u], Zs) 



ujm 
■ -1 



-tg + ao, 



(Ala) 

(Alb) 
(Ale) 

(Aid) 
(Ale) 



From this it is immediately clear that for an individual ray up- 
per turning points coiTespond to the values a = — ^ + 2m, 
while lower turning points are given by a = ^ + 2n7r, where 
n = 0, 1, 2, . . . The horizontal distance between two successive 
bounces on the surface, called a skip-distance, is A(fch,a;) = 
2-Ka = TTCJ m/kf^g = Vpi^ Trm/g, where Vph — ur/kh is the hori- 
zontal phase-speed. The time, tA, it takes the ray to travel between 
the bounces is linked to the distance: A = gt^/Aivm. 

Since > 0, q is increasing when tg is increasing. With uu 
constant on the phase function graph, and kh = k^o = ks cos 6 
constant along the ray, we need to find kz in ray-coordinates to 
complete the solution dAlalAlcl ) to phase space. This is found from 
the condition that Hamiltonian j2.8a| ) is zero. Indeed from ( 14. 3t and 
dAlbb it follows that k^ — khb cos a/ z, with the choice of the 
square root branch taken care of by Some care should be taken in 
choosing the correct branch of the square root, which is addressed 
via a choice of the inverse sine branch in the definition of qq dAlet 
as follows. For given Zs , co and 6, there will be two solutions for pa- 
rameter qo(6»): q+ G [-f ,f] andQo = 7r-a+ G [|,^].The 
former clearly corresponds to rays going down from the source, z 
increasing with tg, kzo > 0, and the latter to rays going up (z de- 
creasing, kzo < 0). Thus, is chosen when 9 > 0, and when 
9 < 0. To ensure smooth dependence of ao on 9, one has to join 
the two solutions at ^ = 0, which corresponds to kzo = 0, i.e. the 
ray that has source depth, Zs, as one of its turning points. In view 
of the Proposition [T] and Corollary [T]in Section 1431 where a parti- 
tion depth is deduced and determined for polytrope model, depend- 
ing on whether the source is located above or below the partition 
depths ze, this turning point will be either upper (qq = f + 2nn) 
or lower (ao = — f + 2nn). 



To treat this mathematically, let us define 



The condition that waves are not evanescent, fef > 



UJ > 



ijJac{zs), implies that ^ G (1,00). When Zs 
helpful to rewrite 



Ze, ^ = 2. It is 



2 cos2 9 



l^((e-2)V4(e-l)sin^e) 



4 cos* 61 

Note that b = only when both ^ = 2 and 61 = 0. When 6/0, 
6 sin ao — Zs — a and b cos ao = Zskzo /kh ~ Zs tan 9, so that 



tan ao — 



1 



|tan( 



cot 6* 



tan 9 2 sin 261 ' 

Consider the behaviour at S = 0, as well as 9 approaching edge 
points, — f , f , that can be potentially problematic. There are three 
cases: 



• Zs> ZE,^e (1,2) 



6* = corresponds to lower turning 



corresponds to upper turning 



. Then, ao 



as 9 



point. Thus chose ao(^ = 0) = ^. Then, ao 

and ao '^'^ as 61 -> + ; 

• Zs < ZE,^ ^ (2,00) => 

point. Thus chose ao {9 — 0) 
f + ;andao -f" us 9 ^ -f + ; 

• Zs — ze,£, = 2,tanao = — tan 6, when ^ 7^ 0. For ^ = 
the value of ao is not defined, but as a = z_b, let us choose 
ao(6' = 0) = with ao -f + as 61 ^ f "; and ao ^ f + as 



Therefore, using formulas for 6,tanao and cosao, the derivative 
of ao by 9 can be expressed in terms of ^: 

dao _ 2(^ - 2) - 4(5 - 1) sin^ 9 
~d9 



= -1 



(C- 2)2 +4(5-1) sin2( 
C(^-2) 



(5-2)2+4(C-l)sin^ 



(A3) 



In the case when Zs = ze,^ = 2, the above formula takes form 
= — 1. For fixed frequency and source depth, taking into ac- 
count the information above, relationship ( IA3t can be integrated to 
find 



ao(9) = I 



6* + tan ^ ( ^ tan 61 



9 + tan 



U-2 



tan ( 



Zs >ZE,i£ (1,2), 
Zs < ZE,i£ (2,cxj) 

Zs = ze,£, ^2. 

(A4) 

Representative plots for ao as function of 9 are shown in Figure|7] 
A2 Caustics 

Equations JAIalAIdt together with with 



kh ~ ksiuj, Zs) cos 9, 
b cos a 

rZz = rZh , 
Z 

Lo — const, 



(A5a) 
(A5b) 
(A5c) 



define the three-dimensional graph of ip lying in Q,M in terms of 
ray coordinates [9, to,tg). As t — tg + to and since neither x nor 
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z depend explicitly on to, the Jacobian of the coordinate transfonn 
can be written as 



V 



d{x, z, t) 
d{e,to,tg) 



-kh 



— 

^ n —k 

ae u u 

1 1 



d{x, z) 

'WJ7) 



(A6) 



Note that D = when the wave- vector {dx/dtg, dz/dtg) is tan- 
gent to the wavefront, i.e. collinear to [dx/dO, dz/dO) . 

When 6 £ (— ^, -l) fej, / 0, a can be used instead of tg 
as a parameter along the ray: 

d[x,z,t) d{e,to,a) ^ khgd[x,z) ^^^^ 
d{9,to,a) d{9,to,tg) Lom d{6,a)' 



This is then integrated to reconstruct the phase function along the 
ray 

m(m + 2) da' 



Lp= ipo ■ 



4fch 
1 



\Jm(m + 2) 



a + b sin a 
tan^ 



a a b 
tan — + ——= 

ZuZl Z \/ ZuZl 



tan 



a ao , b 

—= tan — + ——= 

'ZuZl Z 'w ZuZl 



(All) 



This paper has been typeset from a TgX/ file prepared by the 
author. 



Let us consider the general case z ^ ze- Making use of the 
formulas 



2a tan t 



— ^-fc.tane, - 

db taxiO , 2 ,2x , 

— = — ^'"^ ^ when Zs ze, 

da \. a , 



the partial derivatives with respect to take-off angle are evaluated: 

dx 



dz 
06 



-{a — ao) tanS(a — 6 sin a) 
+ 2s(l + tan^ t 
+ osma-— — 



tan 6* 2 , ,2x 
— - — (a +0 ) cos Of 



^ tan S , 2 ,2\ • 
-Za tan 6* H — (a + b ) sm a 



(a — ao)&cosatan6' + 6 cos a 



dap 
dd 



dx _kh£ 

dtg ujm 

dz _khg 
dta 



(a + 6 sin a), 



= bcosa 



(A9a) 

(A9b) 
(A9c) 
(A9d) 



Then using a as a parameter along the ray define the scaled jaco- 
bian 



J(9, a) — - — V — feces QX 
khg 

( 3a 
X I 2a(a — ao) tanS — a-— - + 2:5(1 + tan^ ( 



tan Q 



{((I + 3fe^a) sin a + ia'b + }? 



(A 10) 



A3 Reconstructing phase function 

Using a as a parameter along the ray and taking into account ( I2.8al . 
phase function derivative along the ray ( 12.10) takes the form 

d</3 ^ dx ^ dz dt z cj^j, 

da ^ da " da da kh 



